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ABSTRACT 

We systematically study the effects of varying the starting redshift Zi for cosmological simulations 
in the highly non-linear regime. Our primary focus lies with the (individual) properties of dark matter 
halos - namely the mass, spin, triaxiality, and concentration - where we find that even substantial 
variations in Zi leave only a small imprint, at least for the probed mass range M G [10^°, 10^^]/i~^ Mq 
and when investigated at redshift z = 0. We further compare simulations started by using the 
standard Zel'dovich approximation to runs based upon initial conditions produced with second order 
Lagrangian perturbation theory. Here we observe the same phenomenon, i.e. that differences in the 
studied (internal) properties of dark matter haloes are practically undetectable. These findings are 
(for the probed mass range) in agreement with other work in the literature. We therefore conclude 
that the commonly used technique for setting up cosmological simulations leads to stable results at 
redshift z = for the mass, the spin parameter, the triaxiality, and the concentration of dark matter 
haloes. 

Subject headings: galaxies: halos — cosmology: theory — cosmology: dark matter — methods: n-body 
simulations — methods: numerical 



1. INTRODUCTION 

During the past two decades numerical simulations of 
cosmic structure formation have become a standard tool 
in cosmology. Advances in computational algorithms 
combined with an ever increasing size of the machines 
used to run these simulations have made it possible to 
simulate regions of the universe with unprecedented dy- 
namic range, recently culminating in multi-billion parti- 
cle simulations of ga lactic halos fe.g.. IStadel et al.ll2008t 
ISpringel et al.l [20081 ). However, most of the efforts to- 
wards achieving this goal have been related to refining 
the simulation techniques and studying the differences 
between va rious codes, respectively (cf. Frenk et al. 1999 ; 
Knebeet al. 2000; O'Shca et al. 2005; Rcgaiietal.li2007l: 



Heitmann et al.. ,2007, : ,Agertz et al., .20071 : iTasker et al.l 
20081) . But common to all these simulations - irre- 
spective of the applied code - is the way how the ini- 
tial conditions (ICs) are generated: practically every- 
one in t he field uses the so-called Zel'dovich approx- 
imation (jZel'dovichl Il970l ZA) first employ e d for cos- 
mological simula t ion bylKlvpin fc ShandarinI (|1983D and 
lEfstathio u et al.l ()1985D . despite the notion that other 
methods m ay be mor e suitable a n d accurate, respec- 
tively (e.g. |Penl [l997t ISirkol [2005t iHansen et al] 120071 : 
iJovce et "all I2008D . We like to remind the reader that 
the ZA is a first order Lagrangian perturbation theory 
and hence we will also refer to ZA as "Iptl". There are 
only a few parameters that require specification for this 
method, one of which is the starting redshift Zi] for an 
elaborate discussion of (most of) th e parameters influ - 
encing the ICs we refer the reader to lJovce et al.l ()2008l ). 

There is already a great deal of work out there 
related to the credibility of the ZA w ith respects 
to transients from initial conditions (e.g. IScoccimarrol 
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119981 ICrocce et all 120061: iTatekawa fc Mizunol I2007D. 
the m emory of initial conditions (jCrocce fc Scoccimarrol 
errors in real-space statistical properties o f 
the ZA (e.g. iKnebe fc Domfnguez 2003; Sirko '2005'), 
discreteness efects (iBaertschiger fc Svlos Labini 200ll ; 
iJovce fc Marcos! l2007bl lah as well as descriptions of 
how to accurately generate multi-mass simulations (e.g. 
iNavarro et al.lll997t|Prunet et al.ll2008f) . However, to our 
surprise the literature yet lacks a systematic study of the 
right choice for the starting redshift Zi with respects to 
the internal properties of dark matter halos at today's 
time. 

On the one hand, workers in the field argue that 
the validity of the linear theory - upon which the ZA 
is based ~ is enforced by choosing a starting redshift 
in such a way that the resulting variance of the dis- 
crete density field (from now on referred to as ctb, see 
definition below) i s significantly less than unity (e.g., 
ICrocce et all 120061 : IPrunet et all 120081 ). But how much 
smaller than unity exactly? As a rule of thumb peo- 
ple adopted a va lue of ~ 0.1 — -2 (as pointed out 
by, for instance. ICrocce et all 120061 : iLukic et all 120071 : 
IPrunet et al.ll2008[ ).'' B ut what is th i s choi ce based upon? 

On the other hand, 'Luki c et al.l (|2007f l. for instance, 
provide a useful formula for estimating the starting red- 
shift Zi based upon the requirement that the fundamental 
mode in the simulation box is well in the linear regime 
(Eq.(20) in iLukic et al.]|2007f) : and they additionally 
discuss the necessity to further restrict the initial dis- 
placement of particles to a level that the first crossing of 
trajectories happens several expansion factors after start- 
ing th e simulation (see also IValageasll2002l : ICrocce et al.l 
l2006l) . They argue that it is important to allow for a 
sufficient number of expansion factors between the start- 
ing redshift and the highest redshift of physical signifi- 
cance to ensure that artifacts from the initial conditions 
(e.g. a regular grid structure) are lost. This goes along 

Please note that these references do not argue in favour of 
cFg < < 1 ; they rather refer to this choice as "common practice" . 
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wih the tests presented in iReed et al] ()2003D where it 
is argued that the simulation should be evolved for at 
least an expansion factor of ~ 10 before extracting (mass 
function) measurements. Starting too late (or allowing 
for too few expansion before extracting physical infor- 
mation from the simulation) will delay collapse of the 
first halos acting as seeds for further structure forma- 
tion, e.g. the Zeldovich approximation cannot account 
for shell-crossing wherei n mass piles up as it flows to- 
ward s overdensities (c f., Jenkin s et al.ll200ll : iReed et al.l 
I2003t iHeitmann et al.|[2008 ). But what is "too late"? 

However, one should also not start too early to avoid 
numerical round-off errors and shot noise of the particles 
used to sample the primordi al matter density field (e.g. in 
the case of glass ICs) (e.g., ILukic et al.ll2007D . The main 
objective of this paper is to shed light on this issue and 
study the differences in properties of dark matter halos at 
redshift z = when systematically varying the starting 
redshift Zi and consequently ub- 

We though have to acknowledge that this subject 
has in part b een touc hed upon by ot her people (e.g., 
Jenkins et all 1 2001: R eed et all 120031: [H citmann et al 



2006; 
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Crocce et aL.2006HLukic et al.ll2007; .Tinker et al. 
Jovce et al.ll2008[ ). However, the primary focus of 



these studies was solely the mass (function) of dark mat- 
ter halos and numerical effects in generating the ICs, 
respectiv ely. In that regar ds, it has been shown by, for 
instance. I Reed et al] ()2003D that the starting redshift can 
have a substantial influence on the high-redshift mass 
function. However, such effects should have evolved away 
by lower redshifts since the tiny fraction of matter that 
is in halos at high reds hift is soon incorp orated into clus- 
ters or large groups (|Reed et al.ll2003h. This is actu- 
ally confirmed by the findings of iJenkins et al.l ()2001l ) 
who observed that the mass function at low-z is not 
very sensitive to the starting reds hift. A similar result 
was found by [Tinker et al.l (|2008f ) for halos with mass 
M < 10"/i"i Mq: they compared mass functions ob- 
tained with simulations based upon ZA initial conditions 
but variations in the starting redshift from Zi = 60 to 
Zi = 35. However, the difference increased for halos more 
massive than lO^^ft."^ Mq and could become as large as 
10-20%, nevertheless depending on the code and partic- 
ulars of the simulation, respectivel y (cf. their Fig. 14) . 
This goes along with the results of ICrocce et al.l ()2006f ) 
who found a ^ 10% discrepancy at lO^^M© in z = 
mass functions when comparing the standard first-order 
ZA with second-order Lagrangian perturbation theory 
(from now on referred to as "lpt2") for generating ini- 
tial conditions for cosmogical simulations. This may be 
caused by the fact that ZA assumes straight lines for 
particle trajectories whereas lpt2 includes the effects of 
gravitational tides; and the latter are most pronounced 
for regions containing the rarest peaks of largest height 
that ten d to evolve into the largest galaxy clusters at low 
redshift ([Tinker et al. 2008). 

Here we extend and complement all previously men- 
tioned studies by quantifying the impact of the initial 
starting redshift Zi upon the individual properties of grav- 
itationally hound objects at redshift z = other than the 
mass alone. In addition, the same comparison is done 
for simulations started by a 2nd order Lagrangian per- 
turbation met hod (lpt2). This com plements the research 
carried out bv lCrocce et al.l ()2006l ) who presented an in- 



depth investigation of transients from initial conditions 
in cosmological simulations in ZA and lpt2. In contrast 
to that study, we consider scales that are much deeper 
in the non-linear regime. 

2. THE SIMULATIONS 

We ran a series of simulations with the publicly avail- 
able GADGET2 code. All simulations consist of TV = 256^ 
particles in a cubical volume of 25h^^ Mpc side length. 
The cosmology we imposed is compliant with the lat- 
est WMAP results (i.e. » o = 0.28, r?A = 0.7 2, /i = 
0.73,(78 0.76,713 = 0.96; iKomatsu et"alll2008l ). The 
force resolution of our simulations is 2h~^ kpc. As the 
phases in the generation of the ICs were identical across 
all models we are able to cross-compare both individual 
particles as well as individual halos. The latter were iden- 
tified with the MPI-enabled open source halo finder AHF* 
(AMIGA's-Halo-Finder, KnoUmann fc Knebe 2009), w hich 
is based upon the MHF halo finder of I Gill et all (|2004f) : ha- 
los are located as peaks in an adaptively smoothed den- 
sity field of the simulation using an adaptive grid hierar- 
chy based upon a refinement criterion that matches the 
force resolution of the actual simulation (i.e. in our case 
5 particles per cell) ; local potential minima are computed 
for each of these peaks and the set of particles that are 
gravitationally bound to the halo are returned. For ev- 
ery halo we calculate a suite of canonical properties (e.g. 
velocity, mass, spin, shape, concentration, etc.) based 
upon the particles within the virial radius. The virial 
radius i?vir is defined as the point where the density pro- 
file (measured in terms of the cosmological background 
density ph) drops below the virial overdensity Avir, i.e. 
M(< i?vir)/(47ri?3.^/3) = AvirPfc. This threshold Avir is 
based upon the dissipationless spherical top-hat collapse 
model and is a function of both cosmological model and 
time. For the given cosmology it amounts to Avir = 354 
at z = 0. 

We need to mention that we restricted our analysis 
to objects with at least 100 particles corresponding to a 
lower mass cut of M^nin = 7 .2 x 10^h~^ M© an d accord- 
ing to the study presented in lKnollmann fc Kn ebe (200£j) 
we are certain that our halo catalogue is complete at this 
level. The total number of objects found is of order 20000 
with approximately 10 objects with M > lO^^h'^ Mq. 
For a more elaborate discussion of how well our numer- 
ically determined mass functions agree with existing fit- 
ting formulae in the literature w e refer the interested 
reader to the actual AHF paper bv iKnollmann fc Knebd 
(|2009l) . 

The simulations can be split into two distinct sets clas- 
sifying the method applied to displace the particles from 
a regul ar lattice at the initial redshift z^, i.e. using the 
code bv lShkol ((20051) we either apply 

• Iptl: the Zel'dovich approximation or 

• lpt2: 2nd order Lagrangian perturbation theory. 

In each of these sets we varied the starting redshift and 
used the values Zi = 150, 100, 50, and 25 and our conven- 
tion for labelling these runs is summarized in Table [TJ 
That table further lists the starting redshift alongside 



AHF is already freely 

|http : //www . alp . de/Peop le/aknebe 



available 



from 
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TABLE 1 

Simulation labels 
alongside initial 
conditions 
parameters. 



run 




""box 


150-lptl 


150 


0.05 


100-lptl 


100 


0.07 


050-lptl 


50 


0.14 


025-lptl 


25 


0.28 


150-lpt2 


150 


0.05 


100-lpt2 


100 


0.07 


050-lpt2 


50 


0.14 


025-lpt2 


25 


0.28 



the aforementioned rms variance of the matter distribu- 
tion inside the computational volume 

'yl^T^ nk)k^dk (1) 

where fcmin = '2.7: / B represents the fundamental mode 
determined by the box size B and fcmax ttN^/^/B the 
Nyquist frequency that additionally depends on the num- 
ber of particles N used for the initial conditions. 

3. DENSITY FIELD COMPARISON 

Even though the primary focus of this study lies with 
the internal properties of dark matter halos (to be pre- 
sented in Section [4] below) we start with an investigation 
of the effects of the starting redshift (and order of the 
Lagrangian scheme) upon the matter density field. This 
is due to the lack of differences found in Section [4] and we 
hence considered it mandatory to verify our simulations 
against the results obtained by other workers in the field 
who predominantly explored the power spectrum (and 
mass function of halos) in that regards. 

According to the commonly accepted and applied 
method to generate ICs for cosmological simulations we 
refer to model "lptX-050" as our reference model for 
which the rms of the matter field ctb = 0.14 is in the 
commonly used range of 0.1 — 0.2. We therefore com- 
pare all other models in the first set ("Iptl" set) against 
this particular run and in analogy use "lpt2-050" as the 
reference in the lpt2 set. 

3.1. Applied Comparisons 

All models in a given lptl/2 set are compared against 
the reference model started at z = 50. To gauge the 
influence of using first (Iptl/ZA) or second order (lpt2) 
Lagrangian perturbation theory in the generation of the 
ICs we additionally cross-compare runs of these two sets 
against each other that started at the same redshift. In 
summary, this leaves us with ten comparisons summa- 
rized in Tabled] 

3.2. Power Spectra 

We start our comparison with the power spectrum of 
matter density fluctuations. The power spectrum 

P(fc)(x(|<5(fc)|2) , (2) 

where 6{k) is the Fourier transform of the density con- 
trast 5{x), i.e. 



TABLE 2 
Applied Corss-Comparisons. 



comparison 


model #1 


model #2 


#1 


150-lptl 


050-lptl 


#2 


100-lptl 


050-lptl 


#3 


025-lptl 


050-lptl 


#4 


150-lpt2 


050-lpt2 


#5 


100-lpt2 


050-lpt2 


#6 


025-lpt2 


050-lpt2 


#7 


150-lptl 


150-lpt2 


#8 


100-lptl 


100-lpt2 


#9 


050-lptl 


050-lpt2 


#10 


025-lptl 


025-lpt2 



0.10 
0.08 
r 0.06 
2 0.04 



o 0.02 



025-lptl vs. 050-lpt1 




100-lptl vs. 050-lptl 




150-lptl vs. 050-lptl 




025-lpt2 vs. 050-lpt2 




100-lpt2vs. 050-lpt2 




150-lpt2vs. 050-lpt2 




z=0 















k [h/Mpc] 

Fig. 1. — Fractional difference in the power spectrum P{k) re- 
garding the different starting redshifts at redshift z = 0. The 
wavenumber k ranges from the largest mode to the Nyquist fre- 
quency of the particle grid. 




-0.02 • ■ • ■ 1 

1 10 

k [h/Mpc] 

Fig. 2. — Fractional difference in the power spectrum P{k) re- 
garding the different schemes (ZA and lpt2) used to set up the 
initial conditions at redshift z = 0. 



^^^^ " y (2^ ^'^^^^^ ■ ^)'^(^) ' 

is the commonly used statistic to describe the clustering 
of the density field. Here, we compute it by applying 
an FFT on a regular 1 024^ grid using the cloud-in-cell 
(jHocknev fc Eastwoodl fTOSS. CIC) scheme for the mass 
assignment. 

We compare the power spectra of the different runs by 
calculating the fractional difference, e.g. in the case of 
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"150-lptl vs. 050-lptl" the fractional difference in P(k) 
is given by (P(fc)i5o-ipti - ^'(fc)o50-ipti)/-P(A;)o50-ipti- 

The results for the runs started at different initial red- 
shifts Zi with respect to the reference run started at 
Zi = 50 (i.e., comparisons #1 through ^6 in Table [2]) 
are shown in Fig. [TJ The thick and thin lines correspond 
to the ZA and lpt2 runs, respectively. We notice that 
the runs which were set up using the ZA differ more 
strongly (±4%) when changing the respective starting 
redshift than the runs set up with lpt2 (±1%). Or in 
other words, ICs generated using the ZA are more sensi- 
tive to the actual starting redshift (at least with respects 
to the power spectrum analysis) than ICs based upon 
lpt2. We further find, that for the ZA runs, the earlier 
the simulation started the more power we get at 2: = 
compared to the reference model started at z = 50, es- 
pecially at the small-scale/high- /c end. Both of these 
results are compliant with the findings of iCrocce et al.l 
12006) andlMS (|2007f ). However, we like to note that 
given our box size and particle number we probe much 
smaller scales and hence a clustering regime that is highly 
non-linear. 

In Fig. [2] we compare the power spectra of the lpt2 
runs against the power spectra of the Iptl runs for each 
starting redshift Zi at redshift z = (i.e. comparisons 
#7 through #10 of Table We find that, in gen- 
eral, the lpt2 initial conditions lead to more power than 
the ZA initial conditions. Obviously, the effect is big- 
ger the later the simulation started (~ 7% for Zi = 25 
and ^ 1% for Zi = 150), simply because the difference 
between first order and second order Lagrangian pertur- 
bation theory decreases with redshift. And again, this 
observation agrees with the resu lts of other works that 



compared Iptl against lpt2 (e.g. ICrocce et al]l2006l : iMal 
120071: iNishimichi et al.|[2008HHeitmann et al.ll2008f) . 

The bottomline of this sub-section therefore is that - 
whatever we will find in the subsequent study below - our 
simulations reproduce the same trends as found by others 
when varying either the starting redshift or the order 
of the Lagrangian perturbation theory, at least when it 
comes to studying the (power spectrum of the) matter 
density field. 

3.3. Particle Positions 

As all ICs were generated using the same phases we 
are in the advantageous position to compare individual 
particle positions across models. To this extent we apply 
two tests. The first consists of calculating the modulus 
of the difference between those positions 



lArl = 



(4) 



where fi is the position in simulation i and 
ulation j, and the second utilizes the so-called density 
cross-correlation coefficient. 

In Fig. [3] we show for our set of comparisons (cf. Sec- 
tion l3.ip the position difference | Ar| (in units of the force 
resolution, i.e. 2h~^ kpc) for 1% of particles randomly se- 
lected from the total number of particles as a function of 
density as measured at the position of the particle^ in the 

^ The density contrast & = {p — ~j3)/~p has first been calculated by 
assigning the mass of each particle to a regular grid of size 512^. 
Then the grid values have been interpolated back to the particles' 
positions. 
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Fig. 3. — Deviation (normalized to the force resolution) of par- 
ticle coordinates at z = 0. The upper panel presents a random 
sample of 1% of all particles alongside the median Ar in six bins 
in 1 + (5rcf . The lower panel only shows the medians multiplied by 
100, 1, and 0.01, respectively, to avoid crowding. 

reference model started at redshift Zi — 50 (i.e. 050-lptl 
and 050-lpt2, respectively) and hence labelled 5o50-ipti 
in the upper panel and S^ci in the lower panel. The up- 
per panel shows the actual scatter plot for one particular 
comparison (i.e., 050-lptl vs. 150-lptl) together with 
the median while the lower panel shows only the medi- 
ans of Ar in six logarithmically spaced bins. Note that 
the whole particle set has been used to calculate the me- 
dians and to avoid crowding in the figures we multiplied 
the medians by 100, 1, and 0.01, respectively. The er- 
ror bars represent the 25'^ and 75"^ percentiles (slightly 
shifted for each model along the x-axis for clarity). 

We notice the expected trend for |Ar| to increase with 
increasing density contrast, i.e. the differences across 
models are more pronounced in high-density regions. We 
now checked (though not shown here) that the differ- 
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TABLE 3 

Cross-correlation coefficients. 



comparison 


L = 32=* 


L = 64^ 


L = 128^ 


L = 256^ 


L = 512^ 


150-lptl vs. 
100-lptl vs. 
025-lptl vs. 


050-lptl 
050-lptl 
050-lptl 


0.9853 
0.9867 
0.9851 


0.9452 
0.9487 
0.9232 


0.8755 
0.8839 
0.8758 


0.7487 
0.7599 
0.7477 


0.6800 
0.6882 
0.6637 


150-lpt2 vs. 
100-lpt2 vs. 
025-lpt2 vs. 


050-lpt2 
050-lpt2 
050-lpt2 


0.9861 
0.9872 
0.9874 


0.9424 
0.9461 
0.9462 


0.8690 
0.8749 
0.8758 


0.7566 
0.7650 
0.7662 


0.6833 
0.6911 
0.6932 


150-lptl vs. 
100-lptl vs. 
050-lptl vs. 
025-lptl vs. 


150-lpt2 
100-lpt2 
050-lpt2 
025-lpt2 


0.9889 
0.9879 
0.9851 
0.9808 


0.9523 
0.9503 
0.9518 
0.8824 


0.8856 
0.8808 
0.8673 
0.8484 


0.7732 
0.7659 
0.7447 
0.7175 


0.6982 
0.6914 
0.6793 
0.6380 



ences are never larger than the virial radii of the halos 
these particles reside in. However, the observed trend is 
expected: the origin of these deviations is the dynami- 
cal instability of partic le trajectories in the high-de nsity 
regions (e.g. IKnebe et al. 2000; Vall uri et aLll2007D . As 
is well known, the trajectories within virialized systems 
tend to be chaotic and any small differences existing at 
any time moment will tend to grow very fast with time. 
The divergence can thus be expected to be more impor- 
tant in non-linear regions and this explains the trend of 
larger Ar's in denser regions. The differences in the low- 
density regions are substantially smaller, but still larger 
than the the force resolution and hence considered phys- 
ical. However, when investigating underdense regions 
l+(5 < 1 we notice that the median of | Ar| "saturates" at 
approximately the level of 10 x force resolution and hence 
defines the level that marks the minimum expectation for 
the position difference. 

We also observe that the medians do not show con- 
siderable variations when changing the starting redshift. 
Further, the trend for |Ar| to increase with 5 - that ap- 
pears to be independent of Zi - is also of comparable 
strength for Iptl and lpt2. 

However, when comparing Iptl against lpt2 there ap- 
pears to be a drift towards smaller particle position dif- 
ferences (in low-density regions) when moving to higher 
starting redshifts. This is readily explained by the fact 
that at higher redshifts the differences between Iptl and 
lpt2 vanish. Nevertheless, this trend is far less pro- 
nounced in high-density regions. We conclude that the 
differences in P{k) as seen in Fig. [2] therefore stem from 
rather low-density regions. 

We also cross-compared an earlier started model to a 
later started one at its actual starting redshift, e.g. a 
snapshot of run 150-lptl at redshift z = 50 to the ICs of 
simulation 050-lptl. We though chose to not show the 
results as all differences in the positions |Ar| are smaller 
than the force resolution and we hence consider them 
unphysical. 

We like to caution the reader that this particular test of 
investigating the spatial differences | Ar| does not provide 
us with the information which method is better or worse. 
It simply shows that the differences in the final matter 
distribution are practically all at the same level, irre- 
spective of starting redshift and order of the Lagrangian 
perturbation theory used. There is a trend for differences 
in the particles' position to increase in high-density re- 
gions but this trend is the same for cross-comparisons of 



the Iptl, lpt2 and Iptl vs. Ipt2 models, respectively. 

As a second quantitative measure of differences 
we construc t the s o-called cross-correlat i on coefficient 
(IColes et all Il993t ISphnter et all 119981 : iKnebe et all 
[2001 



where i and j specify the different simulations and the 
average is taken over the computational box. 

To compute K, we have calculated the densities 5i,5j 
on a regular mesh using th e triangular-shaped cloud 
(iHo cknev fc Eastwood! 11988'. TSC) density assignment 
scheme, and then used the resulting density field to com- 
pute {5i6j) as well as the corresponding variances ai and 
fjj . We have varied the size of the grid in order to show 
the dependence of the cross-correlation on the smoothing 
scale of the density field. 

The results are shown in Table [3] where it is obvious 
that in all cases the cross-correlation worsens for smaller 
smoothing scales (larger grid sizes). This is compliant 
with the trends seen in Fig. [31 i.e. that |Ar| increases 
with 5: the smaller the smoothing, the smaller struc- 
tures are resolved in the density field. The degraded 
cross-correlation therefore indicates that there are differ- 
ences in locations and/or densities of these small-scale 
structures. If we restrict the correlation analysis to a 
coarse grid, we smooth the particle distribution with a 
fairly large smoothing length and smear out the details 
and differences in the small-scale structure. 

But nevertheless, the absolute values of K are of the 
same order when contrasting different comparisons. One 
may additionally argue that K is slightly higher for the 
lpt2 vs. Ipt2 cross-correlations than for the corresponding 
Iptl vs. Iptl case indicative of the higher order of the lpt2 
scheme, which is in agreement with our findings for the 
power spectra (cf. Fig. [T]). Further - as in the case of 
the power spectra (cf. Fig. [2|) - the cross-correlations for 
Iptl vs. Ipt2 become smaller for earlier starting redshifts, 
e.g. applying a 256"^ grid the K value for 150-lptl vs. 150- 
lpt2 is 8% higher than for 025-lptl vs. 025-lpt2. This is 
exactly the same value as we find for the difference of the 
two power spectra at the corresponding wave number (i.e. 
the Nyquist frequency). 

We also like to note that the value of the density cross- 
correlation coefficient K found here is of the same order 
as for a comparison of the same simulation run with dif- 



6 



7- 1.5 



' ' ' 1 1 1 1 






zl50-lptl 






- 


zlOO-lptl 








z050-lptl 
















zl50-lpL2 










zlOO-lpta 










z050-lpt3 












z025-lpt2 
























< < < 1 1 1 






, , 1 





10 10 10^ 10 

Fig. 4. — Ratio of the mass function at redshift z = of gravita- 
tionally bound objects for all eight models to the respective refer- 
ence model started at 2 = 50. The (Poissonian) error bars measure 
the 3 — cr variance. 

ferent codes unde r matching conditions (cf. Table 2 in 
iKnebe et al. l [200l . 

In summary, we do observe differences in the particle 
positions when changing either the starting redshift or 
the order of the Lagrangian perturbation theory. And 
these differences are obviously more pronounced in high- 
density regions where particle trajectories appear to be 
more "chaotic". However, trends for these differences 
to change with starting redshift are at best marginal, 
irrespective of Iptl or lpt2. 

4. HALO COMPARISON 

The same model comparisons as in the previous sec- 
tion (cf. Table [2]) are now going to be performed using 
the halo catalogues obtained via the halo finder AHF (cf. 
Section [2]). In that regards we will focus o n four quan- 
tities , namely the mass M, the triaxiality (jFranx et al.l 

[TMl) 



T = 



(6) 



where a > b > c are the eigenval ues of the moment of 
inertia tensor, the spin parameter (|Bullock et al.|[200ll ) 



A 



L 



V2 MVR 



(7) 



where L is the absolute value of the angular momentum, 
M the halo mass, R its radius, and V'^ — GM/R, and 
finally the concentration 



R_ 



(8) 



where i?2 measures the point where r'^p{r) peaks, with 
p{r) representing the spherically averaged density profile. 

4.1. Distributions 

We start with inspecting the distribution functions for 
the four quantities under investigation. For the mass 
function n{> M) we decided to plot the ratios with re- 
spects to the reference model started at Zi = 50 whereas 
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Fig. 5. — Probability distribution of the triaxiality parameter 
T (top), the spin parameter A (middle), and the concentration c 
(bottom) for all models. 

for P{T), P{X), and P{c) we chose to plot the actual dis- 
tributions. The results can be viewed in Figs. |4] and [5l 

The dependence of the mass function on the start- 
ing redshift has a lready been investigated by several 
groups before (e.g., Jenkins et al."200r:^ ied et al.ll200a 
iCrocce et a l. 2006; Lukic ct al. 2007; Tink er et al.ll2008D . 

Our own results can be viewed in Fig. [?] where we plot 
the ratio of the mass function of a particular model and 
the fiducial 050-lptl/2 run; the (Poissonian) error bars 
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TABLE 4 

Best-fit parameters 
for the spin 
parameter 
distributions. 



run 


Ao 


O"0 


150-lptl 


0.040 


0.509 


100-lptl 


0.039 


0.501 


050-lptl 


0.040 


0.507 


025-lptl 


0.040 


0.506 


150-lpt2 


0.039 


0.524 


100-lpt2 


0.040 


0.515 


050-lpt2 


0.040 


0.513 


025-lpt2 


0.040 


0.519 



measure the 3-(t variance and are hence proportional to 
3 X iVj^aios where N^l^^^ is the average number of halos 
in the respective bin. 

The (lack of) differences seen in Fig. |4] is consistent 
with the outcome of similar s tudie s in the field (e.g . 
Jenkins et al. 2001; Reed et al. 2003- [Crocce et al.|[2006l: 
Lukic et al..,2007: .Tinker et al...2008i ): most of the pre- 
vious investigations focussed on the very high-mass end 
of the mass function, i.e. 10^'* — 10^^h~^ Mq and found 
transients to be important there. However, that mass 
regime is not probed by our simulations; we are analysing 
halos within a mass range for which none of the previ- 
ously mentioned papers have found any differences either. 

The limited influence of the starting redshift upon 
(low-mass) dark matter halos at redshift z — can also 
be observed for the shape, spin parameter, and concen- 
tration presented in Fig. [5l there are hardly any noti- 
cable differences in the distributions when changing Zi. 
This figure is accompanied by Table |4] for which the spin 
parameter distributions have been fitted to a lognormal 
function 

1 f ln^(A/Ao)\ 

with the two best-fit parameters Ao and (Tq listed in Ta- 
ble m We again note that they are practically indis- 
tinguishable, irrespective of the model and the starting 
redshift Zi. However, we note that the width of the dis- 
tribution as measured by (Tq appears to be marginally 
smaller in the Iptl set. 

Our results from this section indicate that the starting 
redshift has practically no influence on today's attributes 
of dark matter halos - at least not for the properties 
analysed here, namely the mass, the spin parameter, the 
triaxiality, and the concentration and for objects within 
the given mass range 10^° — 10^^h~^ M©. And the same 
holds for the order of the Lagrangian perturbation the- 
ory, i.e. whether Iptl or lpt2 is used to generate the ICs 
has no effect on the particulars of halos at redshift z = 0. 

4.2. Cross- Correlations 

As the simulations were started with identical phases 
we could use the particle IDs to establish a mapping be- 
tween two different simulations. This has been applied 
in Section 13.31 where we presented a direct comparison 
of individual particles (e.g., the spatial difference |Ar|). 
But if we plan to do the same for the halos, we require 



10.0 



a. 



I 




0.1 ^ 

10'° 10" 10'^ 10 
l^oso-ipii [''■"'Mq] 



10.0 



0.1 



_ 150 


-Iptl vs 


050- 


Iptl 




100 


-Iptl vs 


050- 


Iptl 




025 


-Iptl vs 


050- 


Iptl 










■ * # m^' 




_ 150 


-lpt2 vs 


050- 


lpt2 




100 


-lpt2 vs 


050- 


lpt2 




025 


-lpt2 vs 


050- 


lpt2 




i H • u 




I ■ - 


— — 


— 1^ - 


_ 150 


-Iptl vs 


150- 


lpt2 




100 


-Iptl vs 


100- 


lpt2 




050 


-Iptl vs 


050- 


lpt2 




025 


-Iptl vs 


025- 


lpt2 


xO.33 


— inrti^i^ 






a ■ - f ■ . 





10'° 10" 10'^ 10'^ 

Fig. 6. — Cross-correlation of the mass M^ii of individual halos. 
The upper panel presents a random sample of 1% of all particles 
alongside the median mass ratio in six bins in M^^f. We further 
show the curves for mass ratio stemming from a differences in the 
number of particles of ±20 as dashed curves. The lower panel only 
shows the 25*'^ and the 75"^ percentiles as bars around the medians 
multiplied by 3, 1, and 1/3, respectively, to avoid crowding. 

a more sophisticated technique to uniquely cross-identify 
halos amongst different simulations. To this extent we 
utilize a tool that comes with the AHF package and is 
called MergerTree. Originally it serves the purpose of 
identifying corresponding objects in the same simulation 
at different redshifts (and hence the name MergerTree). 
But it can also be applied to simulations of different mod- 
els run with the same initial phases for the initial condi- 
tions like in our case. The MergerTree cross-correlation 
is done by linking objects that share the most common 
particles and has bee n succesfully app lied to similar com- 
parisons before (e.g. iKnebe e t al.^ 2006). 

We consider again mass M, triaxiality T, spin param- 
eter A, and concentration c and the results can be found 
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Fig. 7. — Cross-correlation of the triaxiality parameter T of 
individual hales. The same logic for the panels is used as in Fig.|6] 

in Figs. [BHSl where we plot the ratio of said quantities to 
the reference model against the mass M in that refer- 
ence model. As in Fig. [3] we present in the upper panel 
the actual scatter plot for 1% of the particles alongside 
the median in seven bins while the lower panel this time 
only shows the range from the 25**^ to the 75"^ percentile 
(centered about the median multiplied by 3, 1, and 1/3, 
respectively) marginally shifted for each model on the 
X-axis for clarity. 

We note that (the medians of) the ratios in all in- 
stances are consistent with unity. The "error" bars ~ rep- 
resenting 50% of the halos - are well concentrated around 
unity, in particular in the case of the mass ratio, and do 
not vary strongly with halo mass. We though observe 
that the individual scatter about that median (cf. upper 
panels) increases with decreasing halo mass; ;as a matter 
of fact, the ratio can become substantially large (of order 
two for the mass and order ten for the shape/spin) for 
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Fig. 8. — Cross-correlation of the spin parameter A of individual 
halos. The same logic for the panels is used as in Fig.|6] 

the lowest mass objects. However, we attribute this to 
the differences in the particle positions (in high-density 
regions and hence in and about halos) already noted i n 
Fig. H And as shown in 'KnoU mann h Knebd (|2009f) . 
tiny variations in the particle positions lead to different 
density contours upon which the halo finding algorithm 
of AHF is based. This then entails marginal differences 
in the halo properties that become more apparent at the 
low mass end where halos only consist of few particles. 
To better gauge this explanation we also plot in Fig. [5] 
as dashed lines those curves that mark a difference in 
particle number by 20; and we observe that most of the 
ratios lie between these curves. 

The result of this section is rather remarkable as the 
ratio of properties of cross-identified halos is always con- 
sistent with unity, even if the starting redshift is as small 
as = 25. The observed scatter about unity increases at 
the lower mass end which is naturally explained by vari- 
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Fig. 9. — Cross-correlation of the concentration c of individual 
halos. The same logic for the panels is used as in Fig.|6] 



ations in the number of particles making up the actual 
halo and henc e is rather a peculiarity of the halo find- 
ing algorithm (|Knollmann fc Knebell200"9( ). However, the 
variation of the median as indicated by the error bars in 
the lower panels of Figs. [BHl] clearly shows that the me- 
dians are consistent with unity. 

We are therefore confident that while there are one-to- 
one variations the statistical properties of halos are unaf- 
fected by the starting redshift (cf. Section [471]) . This re- 
sult is again extended to the order of the Lagrangian per- 
turbation theory used for generating the ICs. While we 
still found remnants from transients in the power spec- 
trum (cf. Fig. [2|) and deviations of particle positions (cf. 
Fig. ini) no effects are observed for individual properties 
of dark matter halos at redshift z = anymore. 

5. CONCLUSIONS AND DISCUSSION 



We performed a systematic study of varying the start- 
ing redshift Zi for cosmological simulations. We fur- 
ther used two methods to generate the ICs, namely 
the Zel'dovich approximation (a first order Lagrangian 
perturbation method) and an explicit second order La- 
grangian perturbation code. The resulting snapshots at 
redshift z = were analysed with respects to proper- 
ties of the matter density field as well as the statistical 
and individual properties of dark matter halos. Besides 
of (expected) fluctuations in high-density regions and at 
the low-mass end of the halo population the differences 
are rather marginal. On average, the objects have indis- 
tinguishable properties irrespective of the starting red- 
shift. Surprisingly, even the simulation started as late 
as Zi = 25 for which the rms matter fluctuations are 
as = 0.28 and hence larger than the commonly adopted 
value of < 0.1 — 0.2 gave comparable results to the fidu- 
cial model started at Zi = 50. For the probed mass range 
M e [10^°, 10i3]/i~i Mq, we conclude that (at least at 
low redshift z « 0) the starting redshift has little (if 
any) influence on the (statistical) properties of the halo 
population. Or in other words, the resulting halo prop- 
erties of mass M < 10^^h~^ Mq at redshift z — are 
stable against reasonable variations in the rms matter 
fluctuations as and hence the starting redshift Zi. We 
note that our results are in agreement with the findings 
presented throughout the literature where differences 
have only bee n found o n scales M > 10^'^h~^ Mq (e.g., 
Jenkins et all [2001: Reed et all 120031: iHeitmann et all 



2006"; 'Crocce et al. 20061: iLukic et all 120071: [Tinker et all 

2008; Joyce ct al. 2008). 

We further found that also the order of the Lagrangian 
perturbation theory used to generate the ICs is of little 
relevance. While it leaves an imprint on the matter dis- 
tribution (especially the power spectrum) the study of 
gravitationally bound objects at redshift z = is hardly 
affected, at least when it comes to the spin parameter 
A, the triaxiality T, and the concentration c. This re- 
sult actually goes along with the arguments presented by 
iReed et all (|2003D . i.e. the effects of a "wrong" starting 
redshift should have evolved away by lower redshifts since 
the tiny fraction of matter that is in halos at high redshift 
is soon incor porated into c l uster s or large groups. This 
explai ns whv lJenkins et aP ()2001h as well as lTinker et al.l 
((20081 ) did find that the (low-z) mass function is not very 
sensitive to the starting redshift. 

One concern that may be raised with respects to our 
results is the size of our simulation box and finite vol- 
ume effects, respectively. We concede that the ampli- 
tude of fluctuations at the size of the box is of order 
unity at redshift z = and therefore couplings to modes 
on even larger scales are missing. While this certainly 
affects the halo mass function and the power spectrum 
we though argue that our primary results are not in- 
fluenced. If we were to adjust our mass functions for 
such finite v olume effect s using the recipes outline d in, 
for instance, |R eed ct alj ()2007r i. iLukic et all (2007), or 
iPower fc Knebe (2006) we would need to apply the same 
correction to all our models. But as we only cross- 
compare models amongst each other our findings should 
not be contaminated by such a systematic change. In ad- 
dition, we showed in a previous study that the influence 
of large-scale modes upon the properties of dark matter 
haloes at redshift z = will lead to differences in, for in- 
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stance, the spin parameter of order < 1 5% with the con- 
centr ation being hardly affected at all (| Power fc Knebd 
[2OOI . 

We further acknowledge that our analysis focused on 
simulations primarily analysed at redshift z = 0. And 
our results are not as surprisingly as one may ini- 
tially think as at this time the scales under consid- 
eration are deeply in the non-linear regime and the 
memory of the initial condit ions should have been lost 
(jCrocce fc Scoccimarrol [20061 ). At earlier times the sit- 
uation may be different which explains the results of 
others who found a dependence of the high-mass of 
the multiplicity function especially at high redshift (e.g. 



Reed et all 120031: I Crocce et al.l 120061: iLukic et all 120071: 
Tinker et al. n2oo8i r We leave a more in-depth investi- 
gation of this (using higher resolution simulations) to a 
future study. 
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